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Abstract 

Dual energy computerized tomography has gained great interest because of its abiHty to characterize the chemical 
composition of a material rather than simply providing relative attenuation images as in conventional tomography. The 
purpose of this paper is to introduce a novel polychromatic dual energy processing algorithm with an emphasis on 
detection and characterization of piecewise constant objects embedded in an unknown, cluttered background. Physical 
properties of the objects, specifically the Compton scattering and photoelectric absorption coefficients, are assumed 
to be known with some level of uncertainty. Our approach is based on a level-set representation of the characteristic 
function of the object and encompasses a number of regularization techniques for addressing both the prior information 
we have concerning the physical properties of the object as well as fundamental, physics-based limitations associated 
with our ability to jointly recover the Compton scattering and photoelectric absorption properties of the scene. In 
the absence of an object with appropriate physical properties, our approach returns a null characteristic function and 
thus can be viewed as simultaneously solving the detection and characterization problems. Unlike the vast majority 
of methods which define the level set function non-parametrically, i.e., as a dense set of pixel values), we define 
our level set parametrically via radial basis functions (RBF's) and employ a Gauss-Newton type algorithm for cost 
minimization. Numerical results show that the algorithm successfully detects objects of interest, finds their shape and 
location, and gives a adequate reconstruction of the background. 

Index Terms 

Computed tomography, dual-energy, polychromatic spectrum, parametric level set, inverse problems, iterative 
reconstruction 

I. Introduction 

A conventional computed tomography (CT) imaging system provides a reconstruction of the linear attenuation 
coefficient distribution of an object under investigation. Dual energy techniques, however, allow for more detailed 

O. Semerci and E.L. Miller are with the Department of Electrical and Computer Engineering, Tufts University, Medford, MA, 02155 USA 
e-mail: elmiller@ece.tufts.edu. 



January 18, 2013 



DRAFT 



SUBMITTED TO THE IEEE TRANSACTIONS ON IMAGE PROCESSING, 2011 



2 



chemical characterization of the material using measurements from two distinct X-ray spectra. These methods have 
been applied in a range of application areas including non-destructive material evaluation 1 1 1, medical imaging (21, 
cardiac and coronary imaging O, bone densitometry |4 | and airport/seaport security 0. 

The dual energy idea was initially proposed by Alvarez and Macovsky |6| who modeled the total attenuation 
of X-rays as a linear combination photoelectric absorption and Compton scattering coefficients with corresponding 
empirical energy dependent basis functions. In detail, a cubic polynomial approximation of the polychromatic 
measurement models was used to estimate sinograms of photoelectric and Compton components. Once sinograms 
were estimated, the filtered back projection (FBP) algorithm was applied to obtain photoelectric and Compton 
coefficient reconstructions. The polynomial approximation, however, was found to be accurate only for test subjects 
whose properties were "close enough" to data subjects used to obtain the model coefficients. Also the sinogram 
estimation step amplified the noise leading to erroneous reconstructions especially for the photoelectric component 
161 . Later, basis material decomposition methods, where attenuation coefficients of two or more materials constitute 
the basis set for the total attenuation, were proposed Q, (H. For these methods, the reconstruction problem was 
reduced to estimation of the space-varying weight of each basis material. This approach is accurate when the 
properties in the materials of interest are, in a sense, "spanned" by the spectral characteristics of the basis material 
191 . Optimal selection of the basis set is generally application dependent with a variety of approaches having been 
explored in the literature to date ifTOl , ifTTl . |[T2l . A review of the accuracy of these methods can be found in 
|[T2l . The methods where the estimation of sinograms of basis components is followed by FBP reconstructions of 
the corresponding images are referred as prereconstruction methods. Alternatively in the case of postreconstruction 
type of methods, filtered back projection is used to form separate attenuation images from the high and low energy 
projection data which are then mapped into basis material images ifTSl , Q. It is well-known that both of these 
conventional approaches to the dual energy problem have poor noise properties and are prone to inaccuracies in 
the measurement data as they use FBP flAj . ifTSl , |[T6l . 

In addition to reconstruction methods based on FBP, significant work has been directed toward the use of iterative 
image formation schemes for dual energy CT. Such methods provide the opportunity to more precisely take into 
account the nonlinear relationship between the data and the material properties than is afforded by FBP-based 
algorithms. Moreover, iterative techniques can be tailored to the underlying Poisson nature of the observations and 
provide a natural mechanism for incorporating prior information into the image formation process 1 17 |. Within the 
context of dual energy CT, initial work on iterative methods focused on using the algebraic reconstruction technique 
(ART) and basis material decomposition for iterative beam hardening correction ifTSl , |[T9l . More recently, iterative 
methods based on maximum likelihood and maximum a posteriori (MAP) statistical principles have been developed 
for monochromatic and polychromatic dual energy reconstruction problems lT4l . |[T6l , 1201 arising primarily in the 
medical imaging domain. Such methods have proven to be more accurate than previously proposed iterative and 
FBP based approaches especially in low photon count (i.e., low SNR) scenarios. 

While medical applications have motivated a significant majority of algorithm development for dual energy CT, 
in recent years, there has been growing interest in the application of this imaging modality to problems in airport 



January 18, 2013 



DRAFT 



SUBMITTED TO THE IEEE TRANSACTIONS ON IMAGE PROCESSING, 2011 



3 



and seaport security; specifically the screening of checked luggage O, II2TII . 1221 . 123 1. Unlike the medical imaging 
problem, the range of materials encountered in baggage screening is quite broad. As such, it remains an open 
question as to whether and how material basis-type decompositions of the attenuation coefficient can be applied in 
this context. Indeed, the state-of-the art here is represented by the work of Ying et al. who consider the recovery of 
the Compton and photoelectric coefficients directly in |5|. Following a similar approach to |6|, Ying et al. employ 
a prereconstruction imaging scheme in which the sinograms for the photoelectric and Compton coefficients are 
obtained via the solution of a non-linear constrained optimization problem. As was the case in O, high fidelity 
recovery of the photoelectric coefficient proved challenging due to the domination of the Compton effect relative 
to the photoelectric for this class of problems. Roughly speaking the signal to noise ratio associated with the 
photoelectric component was far smaller than that of the Compton making stable recovery of the former difficult. 

Consideration of these issues associated with the use of dual energy CT for luggage screening provides the 
impetus for the work in this paper. To date, the algorithms developed for the screening problem are rather similar 
to those developed for medical applications despite the fact that the fundamental objectives of these two problems 
are somewhat different. In the medical case, the processing goal is the creation of high resolution images capable 
of supporting accurate diagnoses. For the security problem however, while one may still desire a good "picture," 
an additional goal is the determination as to whether a given object contains illicit materials and if so, how are 
they distributed. Such materials of interest will be embedded in an unknown and generally quite inhomogeneous 
background; however, prior information exists concerning the types of materials that may be of interest. Our approach 
to bringing this and related prior information to bear on the image formation process has not, to the best of our 
knowledge, been considered. 

Motivated by these observation, in this paper we propose a variational approach to image formation adapted to 
the problem of detecting and characterizing "anomalies" in multiple physical parameters located in a cluttered and 
unknown environment. The contributions of this work are as follows. First we propose a new, hybrid model and 
associated inversion methods for the determination of the geometry and contrast associated with these anomalies 
along with a suitable pixel-like representation of the unknown background. Should an anomaly not be present, 
images of the physical properties are provided. In the case where a anomalies do exist, in addition to the images, 
we obtain an explicit representation of the shape, location, and contrast of the objects. Thus, the method we propose 
here simultaneously addresses both the detection and characterization problems. 

Second, we develop a new regularization scheme for problems in which there is a mismatch in the sensitivity of the 
data to the various properties being imaged. As noted above, for the dual energy CT problem, the Compton scatter 
effects tends to dominate the photoelectric effect. Such challenges however are not limited to CT. For example, 
when using diffuse optical tomography for medical imaging problems, we have found sensitivity mismatch when 
attempting to recover multiple chromophore concentration profiles |24|. More generally, the ideas considered here 
are of potential use to a much broader range of so-called joint inversion problems where a collection of heterogeneous 
sensors (acoustic, electromagnetic, optical, mechanical, hydrological) are tasked with developing a single "picture" 
of a given region of space. With each modality sensitive to its own constitutive properties of the medium, the 



January 18, 2013 



DRAFT 



SUBMITTED TO THE IEEE TRANSACTIONS ON IMAGE PROCESSING, 2011 



4 



modeling and regularization ideas considered here represents a step in fusing information in such scenarios. An 
example for such an inversion scheme arises in the geophysical imaging context 1 25 1 where a joint inversion of DC 
resistivity and seismic refraction is performed to model subsurface profiles where cross product of the gradients of 
the reconstructed parameters are used to ensure that the profiles provided by the inversion of each parameter are 
in total agreement. 

While the we present here approach is potentially of broad use, in this paper the specific concern is its application 
to the dual energy CT for luggage screening problem. Here, we view the work as a contribution to the use of iterative 
inversion schemes for this class of problems. Like the MAP approach, image formation is cast as the solution to a 
variational problem comprised of terms reflecting a level of fidelity to the data and prior information. In this context, 
the model we propose and the structure of the prior information are new to the dual energy CT literature. In a bit 
more detail, as was the case in Q, we consider here the recovery of the space-varying structure of the photoelectric 
and Compton scattering parameters We model these quantities in a hybrid manner as the superposition of a 
parametrically defined anomaly (should it exists) and a non-parametric background. 

The anomaly model is comprised of a geometric component that is common in the two images and contrast 
parameters that are specific to the photoelectic and Compton images. In this paper, we employ a newly developed 
parametric level set based representation for the anomalies. Classical level-set methods was developed for modeling 
propagation of curves by Osher and Sethian [1261 and have been widely used in image processing applications 1271 . 
1281 , 1291 as well as for a variety of inverse problems 1301 , ISTl , 1321 , 1331 , 1341 . They provide a topologically 
flexible shape-based formulation, elegantly representing multiple objects with complicated geometries. The number 
of unknowns associated with such a model however is equal to the number of pixels in an image. As such they can 
be difficult to use when solving inverse problems |35|. Recently, level set methods have been defined based on low 
order basis functions expansions for the level set function |36|. In the context of image segmentation and topology 
estimation gradient decent-type methods have been developed for determining the unknown expansion coefficients 
describing the level set function |[37l , |[38l , 1391 , BOl , BTl . An alternative parametric level set model specifically 
adapted to solving inverse problems was considered in |35|. As discussed in |35|, this model is especially attractive 
as there is no need for reinitialization or narrow-banding; all issues frequently encountered with curve evolution 
methods for inverse problems. It suffices to incorporate a regularization term that penalizes anomalies with large 
areas. Moreover, because the order of the model is low (on the order of lO's of unknowns) quasi-Newton methods 
become quite feasible for the estimation of the parameters. 

With these advantages, it is this model that forms the basis for the work here. We apply the parametric level set 
idea to the dual energy CT problem where detection of the size number and location of objects plays a crucial role. 
As such we we demonstrate empirically that the detection and localization of these objects becomes possible even 
if the photoelectric image background reconstruction is not very accurate or a FBP based approach would fail to 

^ While a basis-material decomposition could be employed, as discussed previously, it is not clear at the present time which materials are most 
suitable for the screening problem. Moreover, one contribution of the work here is a new approach to stabilize the recovery of the photoelectric 
coefficient. 
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provide reasonable results. 

In addition to this approach for describing the shapes of the anomalies, the second component of our anomaly 
model is the contrast of these object in both the Compton and photoelectric images. Here we exploit the fact that, 
while the nominal background may be quite variable, for applications such as the one motivating this work, one 
is searching for objects of specific composition. In other words, the physical properties (Compton scattering and 
photoelectric coefficients) of any anomaly are known at least to within some degree of precision. We interpret 
this type of information as providing a set of allowable values for the contrast of the object in the photoelectric - 
Compton coefficient space. By definition then, this set of values cannot be assumed by any pixels not in the object; 
i.e., pixels in the background. These two constraints are enforced explicitly in the variational formulation of the 
image formation problem. 

While a parametric model is used to represent the anomaly, given the lack of structure we assume for the 
background, a pixel-type approach is used to model the nominal spatial structure of the photoelectric and Compton 
images. As noted previously, the physics of dual energy CT are such that in noisy situations one can recover the 
Compton image rather well but the photoelectric image tends to be far less accurate El, O, lfT2ll . In an attempt to 
improve this situation, here we take advantage of the fact that the scene being imaged is the same for both parameters. 
More specifically, because physical objects are represented by discontinuities in the constitutive properties of the 
medium, one would expect that edges in both images should be highly correlated. Hence, we introduce a correlation 
based regularization scheme enforcing structural similarity between the gradients of reconstructed images. 

The remainder of this paper is organized as follows. In Section [ll| we provide the polychromatic CT measurement 
model as well as the representation of total attenuation in terms of photoelectric effect and Compton scatter. In 



Section [III] we describe the shape based model via a parametric level set function and formulate the Compton and 
photoelectric images in terms of homogeneous objects of known properties and unknown background pixels. In 
Section [IV] the numerical strategy for the solution is described. We cast the inverse problem in an optimization 
framework in which the solution is the minimizer of a cost function which consists of a data mismatch term 
and regularization terms. To solve the optimization problem, an alternating minimization algorithm for the shape 
and background parameters is applied with Levenberg-Marquardt algorithm employed in each phase. Numerical 
examples are given are given in Section |V| Section |Vl] consists of concluding remarks and thoughts about future 
goals. Finally in the Appendix, we consider a linearized and simplified version of the background reconstruction 
problem and derive lower bounds for the mean square errors (MSE) in photoelectric absorption and Compton scatter 
coefficient estimation. Showing the the bounds are significantly different, we give the reader the intuition about the 
difficulties that arise in photoelectric absorption coefficient reconstruction and explain the motivation that derived 
us to use the correlation based regularization. 

II. Polychromatic Computerized Tomography Formulation 

Typical X-ray sources used in CT applications generate an energy spectra roughly between 20KeV and 140KeV 
ll42i . In this energy range X-ray attenuation is dominated by Compton scatter and photoelectric absorption 01 each 
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of which can be modeled as a product of an energy and material dependent terms as follows 

E) = c{x, y)fKN{E) + p{x, y)fp{E) (1) 

where ii{x^y^E) is the total attenuation, c(x, x) and p{x^y) are the material dependent Compton scatter and 
photoelectric absorption coefficients respectively. The quantity /kn is the Klein-Nishina cross section for Compton 
scattering which is given as: 



JknW — ^ 



2(1 + a) 1 , , 

-i — - - In 1 + 2a) 



— ln(l + 2a) - , 



(2) 



2a ^ ' (l + 2a)2 

where a = £^/510.95KeV. Lastly, fp approximates the energy dependency of the Photoelectric absorption and given 
as 

h = E-\ (3) 

In polychromatic CT the measured quantities are logarithmic projections obtained at several measurement points 
defined formally as: 

P{ct>,x') = -\u^^^. (4) 
Here x') is a Poisson random variable with mean 

Y{^,x')= I S{E)exp{-fKN{E)A,{^,x') 

^ (5) 
-fp{E)A^{<i,,x'))dE^r{^,x') 

where S{E) is X-ray spectrum, (j) is the measurement angle, x' is tilted version of x axis by (/), r{(j)^x') is the 
background signal caused by scatter and detector noise; Ac{(j)^ x') and x') are the Radon transforms for the x- 

ray path (see Fig. [T]) of each measurement of Compton scatter and photoelectric absorption coefficients respectively 
and given as the following: 

A.(^,.') = /c(.,.)^(.cos. + ysin.-.')ciM,, (6) 

Ap{(j)^x') = J p{x^y)S{x cos -\- y sin — x')dxdy. (7) 

We call c{x^y) the Compton and p{x,y) the photoelectric images. The blank scan, Yq, is assumed to be constant 
for all (^, x') and given as 

Yo= [ S{E)dE. (8) 



We define the reconstruction domain P C as a two dimensional rectangular region such that 

V = {{x,y):0<x<l,,0<y<ly}. (9) 

In order to obtain a discrete representation of ^ and ^ we assume that c(x, y) and p(x, y) are piecewise constant 
on each pixel of a regular grid with collection of Np pixels and define the vectors c = [ci, . . . , cnJ^ and p = 
[Pi^ ' ' ' jPNp]^ to be lexicographically ordered collection of pixels values. Now, ^ and ([vj) can be rewritten in 
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Fig. 1. Geometry of the Radon transform 



a matrix form using system matrix A = {a^j} where a^j represents the length of that segment of ray i passing 
through pixel j. 

In this work we assumed parallel beam projections performed for Nq different angles between and 180 degrees. 
For each angle, x-ray intensities on equally spaced points are recorded to give Nm measurements in total. The 
measurement is repeated for low and high energy source spectra to constitute the vector = [m[, m^] which 
consists of 2Nm elements. Using the discrete grid and the system matrix A, the ith measurement for low and high 
energy levels are written as follows: 

[mL]i 

and 

Similar to the continuous case and Yh are vectors of Poisson random variables with means 

[YLh = [ SL{E)exp{-fKN{E)A,,c 

J (12) 

-^(^)A,*p)d^ + rL,i 

and 

^ (13) 

- fp{E)A,,^)dE^rH^i 

where A^^ is the i^^ row of A, Sl{E) and Sh{E), corresponds to low energy and high energy X-ray spectra 
which are shown in Fig. [2] In our work, these spectra with 1 keV energy bins were obtained using SpekCalc, a free 
of charge software program which calculates x-ray spectra from tungsten anode tubes where specifications such as 
tube voltage and filtration thickness are assigned via a friendly user interface 1431 . ll44ll . B31l . Distinct energy levels 
are obtained by alternating the tube voltage between 80 keV and 140 keV. The total number of photons for Sl{E) 



= - In 



In 



[Y. 



L\i 



Yo,L 
0,H 



Yr 



(10) 



(11) 
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is 1.8 X 10^ and for Sh{E) is 3.6 x 10^. The combination of scatter and detector readout noise contributions are 
denoted by rL,i and rH,i- In practice these contributions are not known a priori since they are related to object 
properties (shape and chemical composition) and scanner specification. However, scatter is a smoothly varying 
signal which can be estimated 1461 . BtI I and often assumed known B6ll or characterized by constant primary-to- 
scatter ratio (481, ll49l - On the other hand, detector readout noise characteristics can be approximated by a Gaussion 
distribution ISOl . Here, we assume r^,* and rn^i are additive white Gaussian (AWGN) and call them background 
noise. For a more advanced treatment of scatter see II51L 

in. The Parametric Level-Set Model and Representation of Compton and Photoelectric 

Images 

Compton scatter and photoelectric absorption coefficients are the main interests of our polychromatic CT problem 
and can be regarded as two different images to be reconstructed. Each image is assumed to be formed by piecewise 
constant objects of interest on an unknown background. For a domain Q C V which represents the support of the 
objects of interest, the characteristic function xi^jV) is defined as 

X{x,y) = < (14) 
[ 0, if{x,y)ev\n. 

Now, the Compton and photoelectric images can be written as 

c{x, y) = x(^, y)ca + [1 - x(^, y)] Cb{x, y) (15) 

and 

p{x, y) = x{x, y)pa + [1 - x{x, y)] Pb{x, y). (16) 

Here Ca and Pa are Compton scatter and photoelectric absorption coefficients of the anomaly (object of interest) 
and assumed to be constant; whereas c^^x^y) and p^^x^y) are low order basis expansion representation of the 
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background images given as 

Cb{x, y) = PiBi{x, y) (17) 

i=l 

and 

Ph{x, = ^ aiBi{x, y) (18) 

where {pi} and {ofi}, i = 1,...,A^5 are the sets of Compton and photoelectric expansion coefficients to be 
determined and {Bi{x^y)} is the set of predefined basis functions. The order of the expansion, N^, is determined 
depending on the desired resolution. Pixels or smoother alternatives as radial basis functions can be chosen depending 
on the desired resolution. This approach provides us with the possibility to reduce the number of unknowns 
significantly and is well suited for applications where a coarser representation of the background is acceptable. 
Low order basis expansion formula proved to be successful for monoenergetic x-ray tomography l33l , ITTl . 

The characteristic function x(^, y) is defined to be the zero level set of an Lipschitz continuous object function 
O :V — >R such that 0{x, y) > in 0{x, y) < in n\D and 0{x, y) = in dn. Using 0{x, y), x(^, y) is 
written as 

x{x,y)=H{0{x,y)) (19) 

where H is the step function. In practice, and its derivative 6^ which are smooth approximations of the step 
function and Dirac delta function respectively, are used 1521 . For e G we have the following: 

1 (l + | + isin(^)) ifN<e 

1 if a; > e (20) 

if X < -e 



and 



^(l + 5m(^)) if|x|<e 



5,{x)={ ' '- (21) 

otherwise. 

This kind of implicit (i.e. level set) representation of regions is a well-known approach successfully applied to image 
processing 1271 as well as image formation problems 1531 . The object function y), is represented parametrically 
using a predefined basis set as 

L 

0{x, y) = ^ aiPi{x, y) (22) 

where a^'s are the weight coefficients whereas Pi{x^y) are the functions which belong to the basis set of V = 
{pi,P2, • • • The shape estimation problem then is reduced to the determination of a set of expansion coefficients 
using a quasi-Newton method (see Section 5). Consequently, we eliminate the requirement of a reinitialization 
process as well as implementation of narrow band methods which are essential for the standard approach where a 
signed distance function is used for the level set function and a curve evolution is performed. Additionally, well 
known properties of the level set approach of such as topological flexibility and capability to represent multiple 
objects is maintained 
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In this work, we used exponential radial basis functions (RBF's) for the background basis set {Bi{x^y)\i 
1, . . . , Ni)} and object function basis set {pi{x^ = 1, . . . , L}. Exponential RBFs are defined as 



^,(r)=exp^- " J. (23) 

Here : 7^^ ^ 7^ is the i^^ RBF, r = {x^y) is the Cartesian coordinates vector, ri is the center of i^^ basis 
function and a is the width which is fixed for all RBFs in the basis set. Once the number of elements in the basis 
set is determined a's can be determined by a least squares fit to simple geometries 



Let us gather all the parameters of our model in vector 0^ = 



where a = [ai, . . . , ai,] , 

OL = [«!,..., q^atJ^, and (3 = . . . , /^atJ^. Now, for the set of model parameters 6, the high and low energy 
projection models given in ( [lo| and |ll| can be expressed in a operator form: K{6) = [K {6)J^ ^ K {6)Jj]^ . We 
define the matrix B G "JZ^p^^^ as the discretized basis matrix where [Bj^j is the value of Bj{x^ y) at the center of 
fi^ pixel. Now, the vector of Compton background and photoelectric background images can be expressed as B/3 
and Ba respectively. Finally, we denote the lexicographical order of discretized object function 0{x^y) with the 
vector O e U^^. 

IV. Reconstruction Algorithm 

In our variational framework, we seek for an estimate 6 of the unknown parameter vector minimizing an 
objective functional F{6) which is composed of a data misfit term and regularization terms which incorporate prior 
information, subject to a set of non-linear constraints as 

F{e,X) = \{K{e)-rafi:{K{e)- m) 

minimize 

subject to (pa,Ca) G r 

([B/3]„[Ba]0 ^r, for z = l,...,iV^. 
Here the constraints account for our prior knowledge about the physical properties of the objects. We interpret 
this type of prior information as providing a subset F of allowable values for the contrast of the object in the 
Ca — Pa parameter space which also corresponds to pairs of values that cannot be found in the reconstruction 
of the background. While in general the set of allowable parameter values could be comprised of a number of 
disconnected regions (as would be the case where there were a number of possible chemical compounds of interest 
each characterized by its own uncertainty set), for simplicity here we consider the case where the constraint set is 
elliptical in shape. We define F G 7^^ as an ellipse with the center (co,po) major axis dp and minor axis (Jc (see Fig.js] 
for a symbolic representation). Now, we can write two sets of non-linear inequality constraints gi{pa-, Ca) '-IZ^ 
and ^) • 'Tl^^^ 1Z^p regarding the object and background contrast values respectively as 

giyPa^Ca) = ^ + 5 1 (25) 
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and 

([B/3],-co)^ ([Ba],-po)' 



(26) 



, for i = 1, . . . , Np. 



The optimization problem with these inequahty constraints can be rewritten as 

F{e,\) = l{K{e)-infi:{K{e)-m) 

minimize 

^R{0,X) ^27) 
subject to gi{pa,Ca) < 0, 
g^{(3,a)<0. 

X-ray interactions, hence the photon count at the detector, is governed by a Poisson process f42]. Statistical 
inversion algorithms, where the exact log-likelihood of Poission distribution is used, have been successfully applied 
to medical imaging problems L2QL 1161 , BSl . Sauer and Bouman 1551 derived the weighted least squares data 
misfit term given in ([24]) and ( [27] ) as a quadratic approximation of the Poisson log likelihood function which is 
successfully applied to dual energy CT problem |[T4ll . In their approach the error terms are weighted by the number 
of counts at the detector such that X) = diag(m). In transmission tomography measurements with low intensities 
(low photon counts) correspond to rays coincide with highly attenuating objects and have low SNR. Therefore, the 
measurements with high intensities are considered to be more reliable and highly weighted. From a different point 
of view, if the photon count is large as will be the case here, the Poisson distribution can be approximated by a 
Gaussian distribution ISOl , 1561 . Furthermore, energy integrating detectors used in commercial CT scanners cause 
the statistics to deviate from Poisson distribution 1571 and their noise characteristics are well approximated by a 
Gaussian distribution |[T4l . 1581 . In that case one can assume that the total measurement noise is additive Gaussian 
with zero mean and invertible co variance matrix Fj^Qj^g and set X) equal to T"^.^^. In this work we assume that 
each measurement is independent and of equal quality; hence set T, = al. This approach had been successfully 
applied to limited data CT 1591 application and we feel it is reasonable for this work since artifacts which are 
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Fig. 4. Comparison of mass attenuation coefficients due to Compton scatter and photo electric effect for water. 



caused by high density materials are complications that we prefer to consider at a later time (See Section [Vl|). 

The regularizer R{6) has two terms which penalizes over- sized objects of interest and allows for the use of 
photoelectric absorption-Compton Scatter basis model. 

R{e,X) = Ri{a)^R2{cx,^). (28) 

First is the commonly used penalty term 1271 . 1311 that encourages objects of interest to have a small area: 

Ri{a) = X,\\H{0)\U. (29) 

As stated in fl) and ([3]), due to the E energy dependency of photoelectric absorption, especially for photon 
energies greater than 50 KeV, the contribution of the photoelectric absorption to the total attenuation is very 
small compared to that of Compton scatter |48|. Fig. |4] demonstrates a comparison of the contributions of both 
phenomenon to the total x-ray attenuation for water as a function of energy. The cross section data for this figure 
is taken from Xcom: Photon cross section database |60|. For instance at 80 KeV Compton scatter cross section 
is 1.7 X 10~^cm^/g whereas photoelectric absorption cross section is 5.77 x 10~^cm^/g (i.e. roughly 30 times 
smaller). Hence, there is severe mismatch between sensitivities of the measurement model to these components. 
Therefore, for noisy data, reconstruction of the photoelectric component requires additional effort especially as the 
total attenuation of the scene increases. In the Appendix, we perform a sensitivity analysis of linearized version of 
the background reconstruction problem and show that the lower bound of the error in the estimation of photoelectric 
absorption coefficient is significantly larger than that of the Compton scatter coefficient. One way to ameliorate such 
a situation is to introduce an appropriate regularization scheme. As we are reconstructing two different parameters 
(Compton scatter and photoelectric absorption coefficients) of the same object, one would expect spatial co- variation 



January 18, 2013 



DRAFT 



SUBMITTED TO THE IEEE TRANSACTIONS ON IMAGE PROCESSING, 2011 



13 



between reconstructed images of these parameters. That is to say orientation of the boundaries should be similar 
although intensity variations within individual images may differ. We exploit this a priori information by enforcing 
similarity between gradient vectors of the images of these parameters. This kind of approach was previously used in 
the context of geophysical imaging where Gallardo et al performed a joint inversion of DC resistivity and seismic 
refraction to characterize 2D subsurface profiles |[25l , I6TII . Their approach to enforcing structural similarity via a 
constraint in the optimization process required that the cross product of the gradient vectors associated with the 
two properties should vanish. 

In our case, however, we seek structural co-variation (i.e., two images should be similar in some sense but might 
not be identical) between reconstructed images. Our aim is to incorporate this a priori information to regularize 
our problem in a variational sense rather than imposing a constraint on the optimization process. As the boundary 
of the object of interest is identical in both images (i.e., xi^^v)) it suffices to deal only with background images. 
Herein we introduce the following correlation type of metric 



i?2(a,/3) = A2 



DB/3||i||DBa||i ^ | ' 



[(DB/3)^(DBa)]^ 



(30) 



where D is the gradient matrix defined as D = [D^^ Dy] with the first difference matrices D^^ and in x and y 
direction respectively. It is easily seen that F2 decreases as the correlation of the gradients increase in negative or 
positive direction and vanishes when they are perfectly correlated or anti-correlated. 

The constants {A^}.^^ 2 associated with each term in R{{0) tunes the trade off between the fidelity of the 
solution to the observed data, incorporation of a priori information and regularization. Optimal choice of multiple 
regularization parameters is a challenging problem. L-curve |62|, generalized cross validation(GCV) |[63ll and 
discrepancy principle |64| are the most frequently used methods for optimal choice of multiple regularization 
parameters. Although systematic selection of parameters is provided with these methods, the problem is not 
straightforward and requires running the optimization code several times. In this work, we tuned the regularization 
parameters Ai and A2 by hand such that the values that give the best results judged by the eye were assigned. 
Empirically speaking, for values bigger then an optimum value of A2, reconstruction quality remains the same while 
number iterations to obtain that quality increases, therefore the process becomes more time consuming. 



The solution to ( |27| ) is obtained using an exact penalty approach where the constrained problem is transformed 
to an unconstrained one so that the transformed cost functional has the same minima as the original one. This is 
achieved by adding a term penalizing infeasible regions to the objective functional, then performing minimization 
using an appropriate unconstrained optimization technique 1651 . Specifically, in this work we used the following 
penalty function 

P{r,Ca,Pa:f3,cx) =r ^+(ca,Pa) + (31) 



where f^{x) = max(/(x), 0) and r is the penalty parameter. As differentiability of the penalty term will be required 
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we use the following smooth approximation | 66 | of max(x, 0) function: 

max(x, 0) ^ ^(V + e + x) 



(32) 



where e G 7^^ is small. 

According to exact penalty theorem f65\ If /j. are the Lagrange multipliers of the original problem, for r > 



max(/i^ : i = 1, . . . , 1 + TVp) the local minima of (27) is also the local minima of the following unconstrained 
problem: 



1 



F,{0,\r) = -{K{0) 



m 



m 



(33) 



As opposed to other methods to solve constraint optimization problems such as quadratic penalty method where 
one needs to solve sequence of sub-problems with increasing penalty parameter, it suffices to find one r that will be 
fixed during the minimization process. Finding the Lagrange multipliers, however, is a hard problem. Nonetheless, 
the value of the penalty parameter r can be determined relatively easier than A^'s as the solution is not prone to 
fail with changes in r and it suffices to, heuristically, set it large enough so that the constraints are satisfied 1661 . 
We found this approach satisfactory for the purpose and the minimization algorithm (Levenberg-Marquardt) of this 
work. For an adaptive penalty parameter update strategy for a S/iQP WT\ based method, we refer the reader to the 
work of Byrd et al. (681. 



The minimization of the cost function in (33) is achieved by Levenberg-Marquardt algorithm 11691 . For this aim 
the cost function is rewritten in terms of an error vector e as 

Fp{0) = e{efe{e). (34) 

The error term has the form e{6)^ = [ef , e^^, 63, ] where each term is associated with the corresponding term 
in the cost functional and given as 

€1 = K{e) - m, (35) 



62(a) = V^iHiO), 



and 



€4 



(36) 
(37) 

(38) 



In order to employ Levenberg-Marquardt algorithm, calculation of the Jacobian matrix J, whose rows are the 
derivatives of e{6) with respect to each element in the unknown parameter vector 6, is required. In detail, nth 
column of J is the derivative of e{6) with respect to nth element in the parameter vector 6. 

de{e) 



S{ca,Pa,a,/3,a} 



(39) 
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These derivatives are easily obtained analytically. The details are omitted here, although the Appendix provides 
some intuition. The solution is obtained by updating 6 each iteration as the following 

where h is the solution of the following linear system at each iteration 

(J^J + /il)h = -J^€ with /i > 0. (41) 

Here I is the identity matrix, /i is the damping parameter affecting the size and direction of h and it is adaptively 
tuned during the iteration process fTOll . 

V. Numerical Examples 

To validate our method we performed reconstructions from simulated data according to the data acquisition models 
given in ([TO]) and ^TT) . The background signals Yl and Yh are simulated as white Gaussian noise. Parallel beam 



measurements are simulated for equally spaced 30 angles between and 180 degrees each with 150 source-detector 
pairs. We compare our results with the method proposed in which is based on estimation of two sinograms 
corresponding to Compton scatter and photoelectric absorption then using a FBP method to reconstruct the final 
images. Their work is also based on reconstructing Compton scatter and photoelectric absorption properties with a 
FBP approach; therefore it gives us the opportunity to compare the performance of our method with the standard 
method that is proposed for luggage screening purposes. A Ram-Lak filter multiplied with a Hamming window 
was used in the FBP inversion. We will call this method DEFBP (dual energy FBP). We also applied thresholds to 
the reconstructed images to reveal regions Xt object of interest defined by the region F (see Fig. [Sj. Specifically 
if reconstructed Compton and photoelectric images are Ic and Ip respectively then Xt t)e calculated with the 
Hadamard product as 

Xt = %c o Xp (42) 



where Xc Xp ^re given as 

[^c(p)]*i 



1 if 



(43) 



\.Xc(p)]ij - c{p)o 
otherwise. 

In the minimization process we followed a coordinate decent type algorithm where object boundary (shape) 
parameters ({a^}) are updated while contrast parameters are kept fixed and then contrast parameters (/3, a, c^, 
Pa) are updated while level set parameters are fixed. Convergence is checked for each individual cycle using the 
stopping condition liTOil 

||x^-x^-i|| <e(l+||x^-i) or fc>W (44) 

where corresponds to vector of shape or contrast parameters at k^^ iteration, /cmax is maximum number of 
iterations and e is a small, positive number which is taken 10~^. Once both cycles are completed the overall 



convergence is checked using the same criteria given in ( [44| ) on the reconstructed photoelectric and Compton image 
vectors. 
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Fig. 5. The value of the cost function F, versus number of iterations for the example shown in Fig.|6] 

Quantitative accuracy of reconstructions are determined by relative -error which is given as 



El2 



II I 111 ' 



(45) 



Here, I is the reconstruction of either the Compton image or the photoelectric image and I is the corresponding 
true image. The accuracy of the boundary reconstruction % compared to the true boundary % is determined using 
the Dice's coefficient. Assuming % and % are sets of points in space where the points are corresponding pixels of 
the object of interest, Dice's coefficient is given as 

2|xnx| 



(46) 



Here = 1 shows maximum similarity (i.e., the sets are identical) and = shows the sets do not have a 
common point (i.e., the sets are disjoint). 

For all examples, N^j = 676 RBF's with = 6Ax where Ax is the length of one pixel were used for 



background representation given in (17) and (18). This corresponds to a RBF for every 4x4 super-pixel and 676 
unknown coefficients for each background image instead of 10^ which would have been the case for a pixel based 
reconstruction. For the level set function we used L = 144 RBF's with ag = lOAx as the basis set. The level set 
function is initialized randomly whereas background parameters /3 and a are set to 8 x 10 ~^ and 80 respectively. 
These values are chosen so that the contrast of initial background images were small. The regularization parameters, 
Ai and A2 were set to 0.1 and 10 respectively; whereas the penalty parameter, r was set to 10^. 

The first set of phantoms that consist of objects with various shapes and constant intensities on a homogeneous 
background have the size of 20x20 cm and discretized into lOOx 100 pixels. We emphasize that the geometry of the 
photoelectric absorption and Compton scatter phantoms are taken to be same although the contrasts are different 
as the images represent distinct physical properties of the medium. The object of interest parameters that define F 
are set as the following: cq = 0.19, ac = 0.05, po = 5000, ap = 500; and the boomerang shaped object on lower 
left is assigned as an object of interest (see first row of Fig. |6]) with = 0.2 and Pa = 5000. These parameters 
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Fig. 6. Simulation with the first phantom with 60dB background noise. First row: Ground truth images. Second row: initial images. Third row: 
Reconstruction after 3 cycles of shape and contrast updates. Fourth row: Final reconstructions after 11 cycles. First column: Compton image. 
Second column: Photoelectric image. Third column: Characteristic function, of the object. The solid yellow is the true object while the thin 
line represents the estimated boundary. 



are in the range that correspond to real materials that may be present in a luggage and are chosen generically to 
test our algorithm on objects with moderate attenuation properties. 

In the first example the first phantom was used and dual enery measurements are simulated for 60dB background 
noise. Fig. [S] shows the value of the cost function Fp{0^ A, r) during minimization process. Fig. |6] shows the actual 
phantoms for Compton and photoelectric images, initializations for both images as well as the object boundary, 
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Fig. 7. Simulation with the first phantom with 60dB background noise. First row: Proposed method without regularizer, R2. Second 
row: DEFPB reconstruction. Third Row: DEFBP reconstruction without background noise. First column: Compton image. Second column: 
Photoelectric image. Third column: Characteristic function, of the object. The solid yellow is the true object while the thin line represents 
the estimated boundary. 

TABLE I 

Error analysis for the simulation using the first phantom with 60dB background noise 



Method 



Compton photoelectric 



Proposed method 0.0741 0.0873 0.9387 

Proposed method without R2 0.0818 0.1977 0.9283 

DEFBP 0.4166 129.57 

DEFBP without b.ground noise. 0.1025 0.3011 0.4229 



reconstructions after 3 cycles of shape and contrast updates and the end results after 11 cycles. In the rightmost 
column, the object boundaries are shown with black solid lines where the yellow object shows the support of the 
true object of interest. Fig. [7] shows the reconstruction using the proposed method without using the correlation 
based regularizer R2 and DEFB reconstructions with 60dB background noise and without background noise. As 
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Fig. 8. Simulation with the first phantom with 40dB background noise. First row: Proposed method. Second row: Proposed method 
without regularizer, R2. Third row: DEFPB reconstruction. First column: Compton image. Second column: Photoelectric image. Third column: 
Characteristic function, of the object. The solid yellow is the true object while the thin line represents the estimated boundary. 

TABLE II 

Error analysis for the simulation using the first phantom with 40dB background noise 











Method 


Compton 


photoelectric 




Proposed method 


0.0908 


0.1014 


0.9455 


Proposed method without R2 


0.0818 


0.1886 


0.9113 


DEFBP 


1.7684 


293.6438 


0.0118 



observed from Fig. |6] and Fig. [7] our method provides an accurate reconstruction of the boundaries of the object of 
interest as well as a reasonable reconstruction of the background images where different objects are clearly visible. 
On the other hand, DEFBP method provides reasonable reconstructions only in the absence of background noise 
especially for photoelectric image. It is also observed that using the correlation regularization significantly reduces 
error in the photoelectric image reconstruction. Quantitative evaluation of the reconstructions in the first example 
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Fig. 9. Simulation with the duffle bag phantom with 40dB background noise. First row: Ground truth images. Second row: Proposed method. 
Third row: Proposed method without regularizer, R2. Third Row: Proposed method without regularizer, R2. First column: Compton image. 
Second column: Photoelectric image. Third column: Characteristic function, of the object. The solid yellow is the true object while the thin 
line represents the estimated boundary. 

is given in Table |l] The correlation regularization reduces for the photoelectric image reconstruction roughly 
by a factor of three. 

For the second example the background noise is increased so that SNR was 40dB. Results are demonstrated in 
Fig. [8] Error in the photoelectric image reconstruction increased compared to the previous example where SNR 
was 60dB. Nevertheless, the object boundaries are reconstructed quite accurately. Similar to previous case error 
levels in the reconstructions obtained with DEFBP method are high. Indeed, even for the Compton scatter image, 
the objects in the scene can hardly be distinguished by eye. See Table |Il| for tabulated error values. 

To construct the second phantom we took a DICOM image obtained from a CT scan of a duffel bag and imposed 
Compton scatter and photoelectric absorption coefficients so that the circular object close to middle was an object 
of interest and background was composed of low density clutter. The object of interest characteristics are set as 
the following: cq = 0.3, dc = 0.05, po = 5000, ap = 500. Other parameters of the reconstruction scheme are kept 
the same as the previous example. We performed the simulation with 40dB background noise. Results are shown 
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Fig. 10. Simulation with the duffle bag phantom with 40dB background noise. First row: DEFPB reconstruction. Second row: DEFPB 
reconstruction without background noise. First column: Compton image. Second column: Photoelectric image. Third column: Characteristic 
function, of the object. The solid yellow is the true object while the thin line represents the estimated boundary. 

TABLE III 

Error analysis for the experiment with the duffle bag phantom and 40dB background noise 



Method 



Compton photoelectric 



Proposed method 0.1365 0.1913 0.7468 

Proposed method without i?2 0.1311 0.2971 0.6325 

DEFBP 1.6702 678.99 

DEFBP without b.ground noise. 0.1099 0.2258 0.5026 



in Fig. [9| Fig. [T0| and Table [m 



In the last example we intend to validate our method in the case where there is no object of interest in the scene. 
We repeated the first simulation keeping all the parameters but the definition of object of interest characteristics 
unchanged. For the definition of F we used the following values: cq = 0.12, ac = 0.05, po = 3000, ap = 500. 



The results are shown in Fig. 11 In this example values for Compton and photoelectric images were 0.1046 



and 0.1498 respectively. As a result of the absence of an object of interest in the scene the algorithm returns a null 



characteristic function. The results are demonstrated in Fig. 11 and Table IV 



VI. Conclusions 

A novel semi-shape based polychromatic dual energy CT algorithm with an emphasis on detection of objects 
whose physical characteristics in terms of Compton scatter and photoelectric absorption coefficients are known 
to some degree is developed. The a priori knowledge about the object characteristics is incorporated into the 
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Fig. 11. Simulation with the first phantom with 60dB background noise in the absence of object of interest. First column: Compton image. 
Second column: Photoelectric image. Third column: Characteristic function, of the object. The solid yellow is the true object while the thin 
line represents the estimated boundary. 

TABLE IV 

Error analysis for the experiment with the absence of object of interest and 60dB background noise 



Method 



El2 

Compton photoelectric 



Proposed method 0.1046 



0.1136 



variational framework via constraints where Levenberg-Marquardt algorithm is used for minimization. Parametric 
level set approach where the level set function is described via low order basis expansion is implemented. This shape 
model is capable of representing various kinds of object geometries as shown in Section |V| Along with object of 
interest boundaries the algorithm provides a reasonable reconstruction of the background images. With the proposed 
hybrid method we aim to increase the detection rate of objects of interest in an unknown cluttered background. 
As demonstrated in Section |V] accuracy of the photoelectric image reconstruction decreases as the background 
noise increases; however, object boundaries can still be recovered accurately. Similar to other level set based 
inverse problem methods, segmentation, hence detection, of the object of interest is simultaneously achieved with 
reconstruction of background images. This property is advantageous, for example, in airport security applications. 
Additionally, the use of the proposed correlation based regularization technique improves the reconstruction quality. 
This regularization technique can be used in inverse problem applications where different parameters, whose 
sensitivity to measurement data is vastly different, of the same scene needs to be simultaneously reconstructed 
and spatial similarity between those parameters is expected. 

In future work, we intend to test the algorithm in the presence of metals or similar highly attenuating objects that 
would increase inconsistencies in the measurement data. Low photon count simulations and limited view tomography 
are also of interest of our research and will be addressed in future work. Another goal is to extend the level set 
framework so that multiple object types can be assigned as object of interest and treated individually. Additionally, 
the basis set for the level-set function can be constructed more adaptively. For instance, the widths, the locations 
and the orientation of each element(blob) in the set can be allowed to change during optimization to achieve greater 
flexibility. Consideration of different types of basis functions such as compactly supported radial basis functions and 
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wavelets in the representation of the background images is also a promising extension of this work. Regularization 
parameters are assigned heuristically in this work. For a practically more sound and effective method an automatic 
determination of regularization parameters is desirable. Finally, in order to make our method applicable to real- sized 
problems, we will consider increasing the resolution and phantom diameter, therefore optimizing the source code 
in terms of computational speed-up. 

Appendix 

Sensitivity Analysis of the Background Reconstruction Problem 

In order to give the reader intuition about the problems in reconstructing the photoelectric absorption, we present 
a sensitivity analysis of the background estimation problem. We consider a linearized problem and additive Gaussian 
noise where we analyse the sensitivity of the reconstruction to small chances in the data. We assume no object 
of interest is present and consider the Jacobian matrices that correspond to background parameters for Compton 
scatter and photoelectric absorption. For sake of simplicity we further assume that the Jacobian matrix is full 
ranl|^ and estimate the errors by the relevant diagonal components of the error covariance matrix of the minimum 
variance unbiased (MVU) estimator. We provide lower bounds for errors associated with photoelectric and Compton 
component estimation and show that the bound on the mean square error (MSE) for the latter is significantly smaller. 



Assume we have the measurement model for the low energy spectra Sl{E) as in (10). If we plug (17) and ([18^ 
into (10) and take the derivatives with respect to Pj and aj, we obtain: 

d[mL]i 



[^l] 

- U{E)[A],,^cx)dE 

and 

d[mL\i 



[A] 2* [B] 



J SL{E)fp{E)exp{-fKN{E)[A]i,Bl3 (48) 



[mL]i 

- fp{E)[A]i,Ba)dE. 

The spectra Sl{E) consists of M discrete energy levels, therefore we detine the matrices and Wq, with elements 

[W^]y = SLiEj)exT>{-fKNiEj)[AUBp 

(49) 

-/j,(E,)[A],*Ba) 

and 

= SLiEj)exp{-fKN{Ej)[A]i,Bp 

(50) 

- U{Ej)[A]i,Ba). 

^The assumption is reasonable. For instance, in the first example given in Section [y| the size of the Jacobian matrix associated with the data 
mismatch term is 5640x1354 and its rank was calculated as 1352. 
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Using an appropriate Newton-Cotes numerical integration rule | 71 1 with a weight vector uj we can rewrite ( [47| ) and 
(|48]) as 

Jc = DiR (51) 



and 



(52) 



where Di = diag(W/5diag(ct?)fKN), D2 = diag(Wcdiag(ct?)fp) and = [uil] - ^[A]i^[B]^j. But since Di and 
D2 are diagonal, hence invertable, we can write 



where D = DsD^^. 

Now, we can perform a sensitivity analysis with the following linear system of equations 

Hf + n = (5m 



(53) 



(54) 



where H = [Jc DJc], f = [6/3 Scx]^ . For simplicity here we assume additive white Gaussian noise n ~ 
A/'(0,Icr^). With (54) we investigate how small changes in the measurement data are reflected to the inversion 
process in the presence of Gaussian noise. If H G R^^^, is full rank, the MVU estimator ll72ll is given as 



(55) 



f = (HTH)-iHT5m 
with the associated error covariance matrix, cov(f — f) = E{(f — f)(f — f)^)} given by 

cov(f-f) = a^iU^U)-^ 

H^H H^H 
H^DH RTd^H 

The block matrixes Ai,A2,A3 and A4 are obtained using the matrix inversion lemma ifTTl . The trace of the 
covariance matrix can be used to determine the quality of the estimator. Therefore, only Ai and A4 be of interest 
and each is given as 



(56) 



Ai = H^H - H^DH(H^D^H)-^H^DH 



(57) 



and 



A4 = H^D^H - H^DH(H^H)-iH^DH. (58) 

For the estimates 5/3 and 5a we have the estimation errors Tr(A^^) and Tr(A4 ^) respectively. Now, the aim 
is to provide lower bounds for these traces. To start, define M = DH and consider the reduced singular value 
decomposition (SVD) of M as M = Vm^m^m of H as H = UEV'^. We can write 

A]"^ = V5]"^(I + KK^)"^5]"^V^ (59) 



January 18, 2013 



DRAFT 



SUBMITTED TO THE IEEE TRANSACTIONS ON IMAGE PROCESSING, 2011 



25 



and 

A4 ^ = Vm5]m'(I + KK^)"'5]^'vI^ (60) 

where K = U'^Um ^ R^""^ . Since V is orthogonal, A^^ and + KK'^)-^5]"^ are similar. Likewise, 

A4 ^ and + K.K.^)~^Y,^ are similar. Hence, the traces of A^^ and A4 ^ can be written in terms of these 

similar matrices as 

= tr + KK'^)5]-^] (61) 

and 

= tr [E^^I + K^K)^^] . (62) 



Let us define = [I + KK'^],^ = [I + K'^K],^. Now, we have 

AT AT 



and 



But, we can write 



2 = 1 2=1 



N N 



<72(DH)n'^^(DH) <t2(DH) 



here cri(DH) is the largest singular value of DH. If we plug ([65]) into (64) we get 

A^min(mi) 
'^?(DH) ■ 

Similarly, for Ejs we have 

A^min(mi) 

Next, we use the following inequality ll73ll for the largest singular value of a matrix H 



where, (imax is msix {[Dja). Now, (66) and (65) becomes 



and 



(63) 



2 = 1 ^ ^ ^ 2=1 * ^ ^ 



f 1 _ EL(nM.'^i(DH)) 

a?(DH) Ef=i(nfc^,'Ti(DH)) 
a2(DH) n^i(DH) 
^ iVn^i(DH) _ TV 



-E« > —27^^. (65) 



> —277^- (66) 



al(H)<[||H|U||HU^/^ (67) 

Using ([67]) we can write 

cTi(DH) < ||DH||i||DH|U < d^^ax (||H||i||H|U) (68) 



A^min(mi) 

"-rf^ax(l|H||i||H|U) ^''^ 



A^min(mi) 

^ ||H|U||H|U- ^^^^ 
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TABLE V 

SENSITIVITY Analysis 



Material 


P 


c 


"max 


Water 


4939.2 


0.1907 


4.3841 X 10^ 


Plexiglass 


3670.1 


0.2157 


4.3841 X 10^ 


Aluminium 


72887.5 


0.4547 


7.7815 X 10^ 



The lower bound of Eo, is <imax tim^s bigger than that of Ef^. Our aim is to compare these lower bounds; hence 
at this point, we can be more concrete. For instance, consider a homogeneous 20 x 20 cm area illuminated from 
angle ^ = by a parallel beam of X-rays using the low energy spectra Sl{E). For this set up, we performed 
calculations corresponding three different scenario where the medium is filled with water, plexiglass and aluminium. 
Photoelectric absorption and Compton scatter coefficients |16 | of these materials and calculated <imax values are 
given m Table [a] For instance, in case the medium is water the value ofd ^ is 4.3841 x 10^ which is approximately 
15 times bigger than than the ratio of the actual photoelectric absorption and Compton scatter coefficients. Also 
note that the difference between the lower bounds increases dramatically as the attenuation characteristics of the 
medium increases as in the case of aluminium. 

The lower bound for the error made while estimating the perturbations in photoelectric absorption coefficients is 
significantly bigger than the error associated with the estimation of perturbations in Compton scatter coefficients. 
From this point of view, we expect a similar behavior in the original (nonlinear) problem especially since Gauss- 
Newton (e.g., Levenberg-Marquardt) type algorithms rely on linearization of the problem at the current estimate 
to find the next estimate at each iteration. We can, therefore, conclude that photoelectric image reconstruction is 
more prone to error compared Compton image reconstruction. In this work we tried to address this issue using the 
regularization term given in ( [3Q| . 
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